#include "flow_do.h"
#include "TH1D.h"
#include <vector>
#include <iostream>
#include "TROOT.h"
#include "TMath.h"
#include "TCanvas.h"
#include "TRandom.h"
#include "TF1.h"
#include "TFile.h"
#include "TChain.h"
#include "TClonesArray.h"
#include "TLorentzVector.h"
#include <TStopwatch.h>

using namespace std; // 




int main(int argc, char **argv)
{

		TROOT root("flow","mc");
		if ( argc<6 ) {
				cout<<"++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
				cout<<"usage :  "<<argv[0]<<" <outputFile> <N_target> <N_EP> <Rndm_seed> <Nevt>"<<endl;
				cout<<"++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
				cout << endl << endl;
				exit(1);
		}

		// CONSTANT
		char *outFile = argv[1];
		int N_ch[ND];
		N_ch[kTargetDetector] = atoi(argv[2]);  
		int N_EPdetctor = atoi(argv[3]); 
		int Rndm_seed = atoi(argv[4]);
		int numberEvents= atoi(argv[5]);

		for(int i= 1; i<ND; i++){
 			N_ch[i] = N_EPdetctor;
		}

		TFile *fout = new TFile(outFile,"recreate");

		cout<<"+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
		cout << "Outfile : "<< outFile << endl;
		cout << "Inputs : "<<endl;
		cout <<"\t nPaarticle for Target   = "<< N_ch[0] << endl;
		cout <<"\t nPaarticle for EP Source= "<< N_ch[1] << endl;
		cout <<"\t nEvt = "<< numberEvents  << endl;
		cout<<"+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;

		// ---------- prepareing histograms --------------
		TH1D *hphidist[ND];
		TH1D *hpsidist[NH];
		TH1D *hpsi_caldist[ND][NH];
		TH1D *hdpsi[ND][NH];
		TH1D *hResol[ND][NH];
		TH1D *hVn[ND][NH];
		TH1D *hVn1evt[ND][NH];
		TH1D *hVnReal[ND][NH];
		TH1D *hphidist_measured[ND];

		for(int ih=0; ih<4; ih++){
			hpsidist[ih] = new TH1D(Form("hpsidist%02d", ih), "" , 1024, 0, 2*pi);
		}
		for (int id=0; id<ND; id++) {	
				hphidist[id] = new TH1D(Form("hphidist%02d",id),"phi distribution over all event", 1024, 0, 2*pi);
				hphidist_measured[id] = new TH1D(Form("hphidist_measured%02d",id),"phi_measured distribution over all event", 1024, 0, 2*pi);
			for(int ih=0; ih<4; ih++){

				hpsi_caldist[id][ih] = new TH1D(Form("hpsi_caldist%02d_%02d",id, ih), "" , 1024, -2*pi, 2*pi);
				hdpsi[id][ih] = new TH1D(Form("hdpsi%02d_%02d", id, ih), "" , 1024, -2*pi, 2*pi);
				hResol[id][ih] = new TH1D(Form("hResol%02d_%02d",id, ih), "" , 1024, -1, 1);
				hVn[id][ih] = new TH1D(Form("hVn%02d_%02d",id, ih), "" , 1024, -3, 3);	
				hVn1evt[id][ih] = new TH1D(Form("hVn1evt%02d_%02d",id, ih), "" , 1024, -3, 3);	
				hVnReal[id][ih] = new TH1D(Form("hVnReal%02d_%02d",id, ih), "" , 1024, -3, 3);				
			}
		}






		//***********************************************************************************************//



		// particle generation prepareing 
		TRandom *myRandom = new TRandom(Rndm_seed);
		TStopwatch timer;
		timer.Start();

		double Qx[ND][NH]; 
		double Qy[ND][NH];
		double psi_measured[ND][NH];
		//*****************************************************************************************//
		// ---------------------  event loop start --------------- //
		//*****************************************************************************************//
		for(int ievt =0; ievt<numberEvents; ievt++){
				if( ievt%100==0) {
						cout << " Event: " << ievt << " " << ievt*100.0/numberEvents <<"% "<< endl;
				}
				// random event plane v0, v1 = 0
				double psi_in[NH] = {0, 0, myRandom->Uniform(0,pi), myRandom->Uniform(0,2*pi/3)};
				for (int ih=0; ih<NH ; ih++){
						hpsidist[ih]->Fill(psi_in[ih]);
				}

				// v_n formular //
				TF1 *dNdphi_bulk = new TF1("dNdphi_bulk",
								"[0]*( 1 + 2*[1]*cos(x-[4]) + 2*[2]*cos(2*(x-[5]))+2*[3]*cos(3*(x-[6])))",0,2*pi);
				dNdphi_bulk->SetParameter(0, input_vn[0]); // underyn const(?)
				dNdphi_bulk->SetParameter(1, input_vn[1]);  //v1  - ignore directed flow
				dNdphi_bulk->SetParameter(2, input_vn[2]);  // v2 - from header file 
				dNdphi_bulk->SetParameter(3, input_vn[3]);  // v3 - from header file
				dNdphi_bulk->SetParameter(4, psi_in[1]); // v1 event plane  - 0
				dNdphi_bulk->SetParameter(5, psi_in[2]); // v2 event plane  - random generated 
				dNdphi_bulk->SetParameter(6, psi_in[3]); // v3 event palne  - random generated


				//generatation particle with random angle phi(which followed dNdphi folmuar)
				vector <double> phis[ND];
				vector <double> phis_measured[ND];
				for(int id=0; id<ND; id++){
					phis[id].clear();
						for(int ih=0; ih<NH; ih++){
								Qx[id][ih] =0;
								Qy[id][ih] =0;
								psi_measured[id][ih] =0;
						}
				}
				// initializing done 


				// ************************  inject bulk particle **************************	
				for (int id =0; id<ND; id++) {
					for( int ip =0; ip< N_ch[id]; ip++) {
						phis[id].push_back( dNdphi_bulk->GetRandom() );
					}
				}//detector loop
				// Generation done. ( 0: target(TPC), 1: source(For EP measure) )	


				// ************************  apply detector effect **************************	
				for (int id=0; id<ND; id++){				
					for (int ip=0; ip<phis[id].size(); ip++){
						double m_phi;
						if(det_resol[id]==0){	m_phi=phis[id][ip];}
						else{m_phi=( int(phis[id][ip] / det_resol[id]) * det_resol[id] + det_resol[id]/2. );}

						phis_measured[id].push_back( m_phi );						
					}
				}
				// Generation done. ( 0: target(TPC), 1: source(For EP measure) )	




				// *********************** event plane calculration & Resolution *************
				for (int id=0;id<ND;id++){
					for (int ip=0;ip<phis_measured[id].size();ip++) {
						hphidist[id]->Fill(phis[id][ip]);
						hphidist_measured[id]->Fill(phis_measured[id][ip]);
					}
				}
				for (int id =0; id<ND; id++) {
						for (int ih=2; ih <NH; ih++) {
								for (int ip=0; ip< phis_measured[id].size() ; ip++){
										Qx[id][ih] += TMath::Cos( ih *( phis_measured[id][ip]));
										Qy[id][ih] += TMath::Sin( ih *( phis_measured[id][ip]));
								}
						}
				};

				for (int id =0; id <ND; id++){
						for(int ih=2; ih <NH; ih++){ 
								psi_measured[id][ih] = TMath::ATan2(Qy[id][ih], Qx[id][ih])/ih; // get EP
								hpsi_caldist[id][ih]->Fill(psi_measured[id][ih]);  // fill EP dist
						}
				}
				// event plane calculration done.

			
				for( int id=0; id<ND; id++) {
					for( int ih=2; ih <NH; ih++) {
								double dpsi = psi_measured[id][ih] - psi_in[ih]; 
								hdpsi[id][ih]->Fill( dpsi ); 
								hResol[id][ih]-> Fill( TMath::Cos( ih * dpsi ));
						}
				}
				// resolution calculation done.



				//************************calcuate vn with above particles & measured psi *8*********
				for (int id =1; id<ND; id++) {  
					for (int ih=2; ih<NH; ih++){
						for(int ip=0 ; ip<phis_measured[id].size(); ip++){
								double dphi = phis_measured[kTargetDetector][ip] - psi_measured[id][ih]; 
								hVn1evt[id][ih]->Fill(TMath::Cos(ih * dphi));
								double real_dphi = phis_measured[id][ip] - psi_in[ih];
								hVnReal[id][ih]->Fill(TMath::Cos(ih* real_dphi ));
						}
					}
				}

				for(int id=0; id<ND; id++){
						for(int ih=2; ih<NH; ih++){
								hVn[id][ih]->Fill( hVn1evt[id][ih]->GetMean() );
								hVn1evt[id][ih]->Clear();
						}
				}
		dNdphi_bulk->Delete();
		} // end of evnetp loop

		// ***********************************************
		// Save Histo ***********************************
		// **********************************************
		for(int id=0; id<ND; id++){
				hphidist[id]->Write();
				hphidist_measured[id]->Write();	
				for( int ih=2; ih < NH ; ih++ ){
						hpsi_caldist[id][ih]->Write();
						hdpsi[id][ih]->Write();
						hResol[id][ih]->Write();
						hVn[id][ih]->Write();
						hVnReal[id][ih]->Write();
						
				}			
		} 
		cout <<"Total number of event processed = "<< numberEvents << endl;
		cout <<"Writing histograms into a root file"<<endl;
		fout->cd();
		fout->Close();
		timer.Print();
		cout <<"Sucessfully finished."<<endl;
}
